1. Technical Field
The present invention relates to image processing and more particularly, to applying rapid numerical approximations of convolutions with filters in image processing.
2. Discussion of the Related Art
Many tasks in computer graphics and image processing involve applying large linear translation-invariant (LTI) filters to images. Common examples include low-pass and high-pass filtering of images, and measuring the image response to various filter banks. Some less obvious tasks that can also be accomplished using large LTI filters are demonstrated in FIG. 1: diagram 10 shows reconstructing images 10B by integrating their gradient field 10A; diagram 20 shows fitting a smooth membrane 20B to interpolate a set of boundary values 20A; and diagram 30 shows scattered data interpolation, with 30A as an input and 30B as an output.
While convolution is the most straightforward way of applying an LTI filter to an image, it comes with a high computational cost: O(n2) operations are required to convolve an n-pixel image with a kernel of comparable size. The Fast Fourier Transform (FFT) offers a more efficient, O(n log n) alternative for periodic domains. Other fast approaches have been proposed for certain special cases. For example, one such ad hoc algorithms describes a multiscale approach, which can approximate a convolution with a large Gaussian kernel in O(n) time at hierarchically subsampled locations.
Besides the straightforward O(n2) implementation of convolution, there are certain scenarios where it can be computed more efficiently. Over periodic domains, every convolution operator can be expressed as a circulant Topelitz matrix, which is diagonalized by the Fourier basis. Thus, convolutions with large kernels over periodic domains may be carried out in O(n log n) time using the Fast Fourier Transform.
Convolution with separable 2D kernels, which may be expressed as the outer product of two 1D kernels, can be sped up by first performing a 1D horizontal convolution, followed by one in the vertical direction (or vice versa). Thus, the cost is O(kn), where k is the length of the 1D kernels. Non-separable kernels may be approximated by separable ones using the SVD decomposition of the kernel B-spline filters of various orders and scales may be evaluated at constant time per pixel using repeated integration or generalized integral images.
Many methods have been developed specifically to approximate convolutions with Gaussian kernels, because of their important role. Of particular relevance to this work is the hierarchical discrete correlation scheme proposed. This multiscale scheme approximates a Gaussian pyramid in O(n) time. Since convolving with a Gaussian band limits the signal, interpolating the coarser levels back to the finer ones provides an approximation to a convolution with a large Gaussian kernel in O(n). However, this approximation is accurate only for Gaussians of certain widths.
The aforementioned ideas culminated in the Laplacian Pyramid, and were later connected with wavelets. More specifically, the Laplacian pyramid may be viewed as a high-density wavelet transform. These ideas are also echoed, where a multiscale bilateral filter pyramid is constructed in O(n log n) time. The Laplacian pyramid, as well as various other wavelet transforms decomposes a signal into its low and high frequencies, which was shown to be useful for a variety of analysis and synthesis tasks. In particular, it is possible to approximate the effect of convolutions with large kernels. However, even though these schemes rely on repeated convolution (typically with small kernels), they are not translation-invariant, i.e., if the input is translated, the resulting analysis and synthesis will not differ only by a translation. This is due to the subsampling operations these schemes use for achieving their fast O(n) performance.
One scenario, which gained considerable importance in the past decade, is the recovery of a signal from its convolved form, e.g., an image from its gradient field. This gives rise to a translation-invariant system of linear equations, such as the Poisson Equation. The corresponding matrices are, again, Topelitz matrices, which can be inverted in O(n log n) time, using FFT over periodic domains. However, there are even faster O(n) solvers for handling specific types of such equations over both periodic and non-periodic domains. The multigrid method and hierarchical basis preconditioning achieve linear performance by operating at multiple scales. A state-of-the-art multigrid solver for large gradient domain problems is described, and a GPU implementation for smaller problems is described. Since the matrices in these linear systems are circulant (or nearly circulant, depending on the domain), their inverse is also a circulant convolution matrix. Hence, in principle, the solution can be obtained (or approximated) by convolving the right-hand side with a kernel, e.g., the Green's function in cases of the infinite Poisson equation. Our approach enables accurately approximating the convolution with such kernels, and hence provides a more efficient, and easy to implement alternative for solving this type of equations.
Gradient domain techniques are also extensively used for seamless cloning and stitching, yielding a similar linear system, but with different (Dirichlet) boundary conditions, typically solved over irregularly shaped domains. In this scenario, the solution often boils down to computing a smooth membrane that interpolates a set of given boundary values describes a dedicated quad-tree based solver for such systems, while avoid solving a linear system altogether by using mean-value interpolation.
This approach can also be useful in more general scattered data interpolation scenarios, when there are many data points to be interpolated or when a dense evaluation of the interpolated function is needed. It would be advantageous therefore to provide an efficient approximation for scattered data interpolation where it is important to use large convolution filters that propagate the information to the entire domain.